source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/fdx1/bmsc/'

bmsc_mex <-
list.files(proj.nm, 'No', include.dirs = T, recursive = T, full.names = T) |>
  set_names(nm = c('No1', 'No2')) |>
  Read10X()

bmsc_mex |> glimpse()

sobj <- bmsc_mex |> CreateSeuratObject(min.cells = 3, min.features = 200)

sobj <- sobj |> PercentageFeatureSet('^mt-', col.name = 'mito_ratio')

sobj |> VlnPlot('mito_ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito_ratio < 10)

sobj <- sobj |>
  quick_process_seurat()

sobj <- sobj |>
  mutate(orig.ident = fct_relevel(orig.ident, 'WT'))

# rds checkpoint ---------
sobj |> write_rds(str_c(proj.nm, 'bmsc.wtko.rds'))

sobj <- read_rds('mission/fdx1/bmsc/bmsc.wtko.rds')

sobj |>
  DimPlot(label = T, label.box = T, repel = T, split.by = 'orig.ident',
          label.size = 2, cols = DiscretePalette(36)) +
  theme_jpub(theme_classic)

publish_pdf('bmsc.split.leiden.umap.pdf', width = 100)

sobj |>
  DimPlot(split.by = 'orig.ident', cols = DiscretePalette(36)) +
  theme_jpub(theme_classic)

publish_pdf('bmsc.split.leiden.nolabel.umap.pdf', width = 100)

leiden_frac_fc <- sobj@meta.data |>
  discov_frac_change(group = orig.ident, subtype = seurat_clusters,
                     var.1 = No1, var.2 = No2)

leiden_frac_fc |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('royalblue', 'red2', 'grey'))

# cell type annotation ---------
## hpca --------
hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(hpca, new_label = 'hpca_main')

sobj |> DimPlot(group.by = 'hpca_main')

## Cell Reports 2019 ---------
cr19 <- list(PreAd='Cebpb', AdP='Cxcl1', MSC=c('Il34','Serpina3n'),
             OsP=c('Alpl','Wif1'), PreOC='Spp1', ProOs=c('Bglap','Col1a1'),
             ProCh=c('Mepe','Dmp1'))

sobj |> DotPlot(cr19, cols = 'RdYlBu', cluster.idents = T)

mrk13 <- sobj |> FindConservedMarkers(ident.1 = 13, grouping.var = 'orig.ident',
                             only.pos = T)

mrk13 |> as_tibble(rownames = 'gene') |>
  filter(str_detect(gene, 'Cd\\d'))

mrk12 <- sobj |> FindConservedMarkers(ident.1 = 12, grouping.var = 'orig.ident',
                                      only.pos = T)

mrk12 |> head()

mrk9 <- sobj |> FindConservedMarkers(ident.1 = 9, grouping.var = 'orig.ident',
                                      only.pos = T) |>
  as_tibble(rownames = 'gene')

mrk9

sobj <- sobj |>
  mutate(cell_type = case_when(seurat_clusters == 12 ~ 'Endothelial cells',
                               seurat_clusters == 13 ~ 'T cells',
                               seurat_clusters == 14 ~ 'Pro-Chondrocytes',
                               seurat_clusters == 9 ~ 'Macrophages',
                               seurat_clusters %in% c(1,3,7,10,16) ~ 'Osteo-lineage',
                               .default = 'Adipo-lineage'))

sobj |> DimPlot(group.by = 'cell_type')

sobj |> DimPlot(group.by = 'cell_type', split.by = 'orig.ident')

## all marker ---------
leiden_mrks <- sobj |>
  FindAllMarkers(only.pos = T, logfc.threshold = 1)

leiden_mrks |>
  as_tibble() |>
  slice_head(n = 2, by = cluster) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj),
         gene = fct_inorder(gene)) |>
  ggplot(aes(cluster, gene, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21) +
  #rotate_x_text() +
  scale_fill_distiller()

# cell type logfc -------
celltype_frac_fc <- sobj@meta.data |>
  discov_frac_change(group = orig.ident, subtype = cell_type,
                     var.1 = No2, var.2 = No1)

celltype_frac_fc |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('royalblue', 'red2', 'grey')) +
  theme_bw() +
  labs(x = 'Cell type', title = 'Cell proportion change in KO vs WT')

finetype_frac_fc <- sobj@meta.data |>
  discov_frac_change(group = orig.ident, subtype = fine_type,
                     var.1 = KO, var.2 = WT)

finetype_frac_fc |>
  ggplot(aes(log2fc_frac, subtype, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('royalblue', 'red2', 'grey')) +
  theme_bw() +
  labs(y = 'Cell type', title = 'Cell proportion change in KO vs WT')

finetype_frac_fc |>
  mutate(fine_type = as.character(subtype)) |>
  left_join(x = sobj, y = _) |>
  DimPlot(group.by = 'type', cols = c('royalblue', 'red2', 'grey')) +
  ggtitle('Cell proportion change in KO vs WT')

# adipo-lineage -------
adipo_allmrk <- sobj |>
  filter(cell_type == 'Adipo-lineage') |>
  FindAllMarkers(logfc.threshold = 1, only.pos = T)

adipo_allmrk |>
  as_tibble() |>
  slice_head(n = 6, by = cluster) |>
  write_csv('adipo-lineage.csv')

# osteo-lineage --------
osteo_allmrk <- sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  FindAllMarkers(logfc.threshold = 1, only.pos = T)

osteo_allmrk |>
  as_tibble() |>
  slice_head(n = 6, by = cluster) |>
  write_csv('osteo-lineage.csv')

sobj |>
  DotPlot(c('Runx2','Ocn'), cols = 'RdYlBu', cluster.idents = T)

sobj <- sobj |>
  mutate(fine_type = case_when(seurat_clusters == 16 ~ 'Runx2+ Osteoblast',
                               seurat_clusters == 3 ~ 'Apoe+ Pre-Osteoclast',
                               seurat_clusters == 7 ~ 'Cdh5+ Endothelial-like',
                               seurat_clusters == 10 ~ 'Kif4+ Osteo-Progenitor',
                               seurat_clusters == 1 ~ 'Bglap2+ Pre-Osteoblast',
                               .default = cell_type) |>
           fct_relevel('Kif4+ Osteo-Progenitor','Bglap2+ Pre-Osteoblast',
                       'Runx2+ Osteoblast','Apoe+ Pre-Osteoclast',
                       'Cdh5+ Endothelial-like'))

sobj |>
  DimPlot(group.by = 'fine_type', cols = 'Paired')

sobj |>
  DimPlot(group.by = 'fine_type', split.by = 'orig.ident', cols = 'Paired')

# mTOR/ampk/akt pathway score -------
search_go_term('TOR')

tor_signal <- map_go_gene('GO:0031929', org = 'mouse')

search_go_term('AMPK')

ampk_signal <- map_go_gene('GO:0004679', org = 'mouse')

akt_signal <- read_csv('mission/fdx1/bmsc/akt.chatglm.txt')

sobj <- sobj |>
  AddModuleScore(features = list(tor = tor_signal$SYMBOL,
                                 ampk = ampk_signal$SYMBOL,
                                 akt = akt_signal$gene |> str_to_title()))

sobj <- sobj |>
  mutate(TOR_score = Cluster1, AMPK_score = Cluster2, AKT_score = Cluster3,
         .keep = 'unused')

## mTOR ---------
sobj |>
  DotPlot2d('TOR_score', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  DotPlot2d('TOR_score', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

## AMPK --------
sobj |>
  DotPlot2d('AMPK_score', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  DotPlot2d('AMPK_score', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

## AKT ---------
sobj |>
  DotPlot2d('AKT_score', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  DotPlot2d('AKT_score', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

# cc.genes ---------
cc_signal <- map(cc.genes.updated.2019, str_to_title) |>
  list_c()

sobj <- sobj |> AddModuleScore(features = list(cc_signal), name = 'Cell_cycle') |>
  mutate(Cell_cycle_score = Cell_cycle1, .keep = 'unused')

sobj |>
  DotPlot2d('Cell_cycle_score', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  DotPlot2d('Cell_cycle_score', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  filter(cell_type == 'Osteo-lineage', fine_type != 'Runx2+ Osteoblast') |>
  DotPlot2d('Ccnd1', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

# osteoblast TF ---------
sobj |>
  DotPlot2d('Sp7', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  DotPlot2d('Sp7', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  DotPlot2d('Runx2', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  DotPlot2d('Runx2', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  FindMarkersAcrossVar(split.by = 'fine_type', group.by = 'orig.ident',
                       ident.1 = 'KO', features = c('Sp7','Runx2'))

# Prrx1+ cell ---------
sobj |>
  FeaturePlot('Prrx1')

sobj |>
  FeaturePlot(c('Prrx1','Slc31a1'), split.by = 'orig.ident', order = T) &
  theme_jpub(theme_classic) & NoLegend()

publish_pdf('Prrx1_Slc31a1_WT_KO_umap.pdf', width = 100, height = 100)

sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  write_rds('mission/fdx1/bmsc/bmsc.osteo.wtko.rds')

sobj |>
  DotPlot2d('Prrx1', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  DotPlot2d('Prrx1', group.x = orig.ident, group.y = fine_type) +
  labs(x = 'Sample', y = 'Cell type')

osteo_fine_kovwt <- sobj |>
  filter(cell_type == 'Osteo-lineage') |>
  FindMarkersAcrossVar(split.by = 'fine_type', group.by = 'orig.ident',
                       ident.1 = 'KO')

osteo_fine_kovwt |>
  write_source_csv('prrx1_subtype_KOvsWT_deg')

## GO GSEA --------
library(clusterProfiler)

osteo_fine_kovwt <- read_csv('mission/fdx1/bmsc/results/prrx1_subtype_KOvsWT_deg.csv')

osteo_fine_kovwt |>
  filter(p_val_adj < .05) |>
  dplyr::count(cluster)

osteo_fine_kovwt |>
  filter(p_val_adj < .05, abs(avg_log2FC) > 1) |>
  ggplot(aes(avg_log2FC, cluster)) +
  geom_jitter(height = .1)

osteo_fine_kovwt <- osteo_fine_kovwt |>
  mutate(p_val_adj = ifelse(str_detect(cluster, 'Runx') & p_val < .01,
                            .04, p_val_adj))

osteo_fine_gsego <- osteo_fine_kovwt |>
  batch_enrich_path('Mm')

osteo_fine_orago <- osteo_fine_kovwt |>
  batch_enrich_path('Mm', method = 'ORA')

osteo_fine_gsekegg <- osteo_fine_kovwt |>
  batch_gsea('Mm', path = 'KEGG')

osteo_fine_orakegg <- osteo_fine_kovwt |>
  batch_enrich_path('Mm', path = 'KEGG', method = 'ORA')

osteo_fine_gsego |>
  imap(\(x,y){x@result |> as_tibble() |> mutate(cluster = y)})

osteo_fine_gsekegg |>
  imap(\(x,y){x@result |> as_tibble() |> mutate(cluster = y)})

osteo_fine_orago |>
  imap(\(x,y){x@result |> as_tibble() |> mutate(cluster = y)})

osteo_fine_gsego |>
  iwalk(\(x,y){x@result |>
      write_source_csv(str_glue('{y}_KOvsWT_GSEA_GO'))})

osteo_fine_gsego$`Bglap2+ Pre-Osteoblast`@result |>
  as_tibble() |>
  plot_enrichment(metric = NES, base_col = 'Blues') +
  labs(x = 'NES', title = 'GO GSEA downregulated pathway in Bglap2+ Pre-Osteoblast') +
  theme_jpub

osteo_fine_gsego$`Kif4+ Osteo-Progenitor`@result |>
  as_tibble() |>
  plot_enrichment(metric = NES, base_col = 'Blues') +
  labs(x = 'NES', title = 'GO GSEA downregulated pathway in Kif4+ Osteo-Progenitor') +
  theme_jpub

osteo_fine_gsego$`Runx2+ Osteoblast`@result |>
  as_tibble() |>
  plot_enrichment(metric = NES, base_col = 'Blues') +
  labs(x = 'NES', title = 'GO GSEA downregulated pathway in Runx2+ Osteoblast') +
  theme_jpub

bglap_gsego <- osteo_fine_kovwt |>
  filter(str_detect(cluster, 'Bglap'), p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = str_glue('org.{org}.eg.db'),
        keyType = 'SYMBOL', eps = 0)

apoe_gsego <- osteo_fine_kovwt |>
  filter(str_detect(cluster, 'Apoe'), p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = str_glue('org.{org}.eg.db'),
        keyType = 'SYMBOL', eps = 0)

apoe_gsego@result |>
  as_tibble() |>
  #filter(NES < 0) |>
  plot_enrichment(metric = NES, base_col = 'Blues') +
  labs(x = 'NES', title = 'GO GSEA downregulated pathway in Apoe+ Pre-Osteoclast') +
  theme_jpub

bglap_orago <- osteo_fine_kovwt |>
  filter(str_detect(cluster, 'Bglap'), p_val_adj < .05, avg_log2FC > 1) |>
  pull(gene) |>
  enrichGO(ont = 'ALL', OrgDb = str_glue('org.{org}.eg.db'),
        keyType = 'SYMBOL', readable = T)

bglap_orakegg <-
osteo_fine_kovwt$gene |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID', OrgDb = 'org.Mm.eg.db') |>
  mutate(gene = SYMBOL) |>
  left_join(osteo_fine_kovwt) |>
  as_tibble() |>
  filter(str_detect(cluster, 'Bglap'), p_val_adj < .05, avg_log2FC > 1) |>
  pull(ENTREZID) |>
  head()
  enrichKEGG(organism = ifelse(org == 'Hs', 'hsa', 'mmu'))

sobj |> dplyr::count(fine_type, orig.ident) |>
  filter(str_detect(fine_type, '\\+')) |>
  ggplot(aes(orig.ident, n, fill = orig.ident)) + geom_col() + theme_bw() + 
  labs(title = 'Osteo-lineage cell number', x = 'Sample', fill = 'Sample') +
  geom_text(aes(label = n), nudge_y = 2) +
  facet_wrap(~fine_type)
